import subprocess
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
from scipy import optimize
from IPython.display import Image, display
Генерация входных данных. Функция принимает на вход длину связи СС, угол связи НСС и торсионный угол СС
def input_gen(l,a,t):
inp = f'''!HF RHF 6-31G
* int 0 1
C 0 0 0 0 0 0
C 1 0 0 {l} 0 0
H 1 2 0 1.08439 {a} 0
H 1 2 3 1.08439 111.200 120
H 1 2 3 1.08439 111.200 -120
H 2 1 3 1.08439 111.200 {t}
H 2 1 3 1.08439 111.200 {t-120}
H 2 1 3 1.08439 111.200 {t-240}
*
'''
return inp
Функция для рассчета энергии молекулы исходя из ее z-матрицы
def run_orca(inp):
with open('orca.inp', 'w') as outfile:
outfile.write(inp)
with subprocess.Popen("/srv/databases/orca/orca orca.inp", shell=True, stdout=subprocess.PIPE, stderr=subprocess.PIPE) as p:
out = p.communicate()[0].decode()
for l in out.splitlines():
if 'FINAL SINGLE POINT' in l:
return float(l.strip().split()[4])
return out
Изменяем длину связи от 1.35 до 1.75
x_o_l = np.arange(1.35,1.75,0.02)
y_o_l = [run_orca(input_gen(l,111.200,180)) for l in x_o_l]
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Целевая функция f(x)=k(b-x)^2 + a
errfunc = lambda p, x, y: fitfunc(p, x) - y # Функция ошибки
p0 = [1,1, -79] # Начальные параметры
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_l, y_o_l)) # Функция оптимизации, получаем из нее подобранные параметры
print("Optimized params:", p1)
plt.plot(x_o_l, y_o_l, "ro", x_o_l,fitfunc(p1,x_o_l),"r-",c='blue',alpha=0.5)
plt.xlim(1.2,1.8)
plt.show()
Был хорошо найден минимум энергии, однако зависимость в целом описана неточно. Возможно следует использовать потенциал Морзе.
Изменяем угол НСС от 109 до 113
x_o_a=np.arange(109.0,113.0,0.2)
y_o_a=[run_orca(input_gen(1.52986,a,180)) for a in x_o_a]
fitfunc = lambda p, x: p[0]*pow(p[1]-x,2) + p[2] # Целевая функция по сути такая же, как и для длины связи СС
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1,1, -79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_a, y_o_a))
print("Optimized params:", p1)
plt.plot(x_o_a, y_o_a, "ro", x_o_a, fitfunc(p1,x_o_a), "r-", c='blue', alpha=0.5)
plt.xlim(108,114)
plt.show()
В этом случае оптимизация прошла успешнее
Значения, полученные в результате выполнения данной работы, отличаются от таковых в статье "Development and Testing of a General Amber Force Field". В этой работе использовались силовые поля GAFF - модификация Amber, заточенная под малые органические молекулы. Скорее всего, мы определили значения параметров не очень точно, ибо использовали иные подходы к их изначальной оценке (при вычислении энергии программой Orca). Вк тому же, силовые константы вычисляются эмпирически, поэтому отличающиеся результаты оправданы.
x_o_t=np.arange(-180,180,12)
y_o_t=[run_orca(input_gen(1.52986,111.200,t)) for t in x_o_t]
fitfunc = lambda p, x: p[0]*(1 + np.cos(p[1]*x - p[2]))+p[3] # Тут целевая функция задается как f(x)=k*(1 + cos(a*x - b)) + c
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [0.25,10,0,-79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_t, y_o_t))
print("Optimized params:", p1)
plt.plot(x_o_t, y_o_t, "ro", x_o_t, fitfunc(p1,x_o_t),"r-", c='blue', alpha=0.5)
plt.xlim(-200,200)
plt.show()
Наблюдаются три минимума, соответствующие незаслоненным конформациям
Рассчет зависимости энергии от длины связи в интервале от 1 до 3
x_o_l_2=np.arange(1,3,0.1)
y_o_l_2=[run_orca(input_gen(l,111.200,180)) for l in x_o_l_2]
fitfunc = lambda p, x: p[0]*((1-np.exp(-p[1]*(x-p[2])))**2)+p[3] # На этот раз описываем потенциалом Морзе, f(x)=k*(1 - exp(-a(x - b)))^2 + c
errfunc = lambda p, x, y: fitfunc(p, x) - y
p0 = [1,1,1,-79]
p1, success = optimize.leastsq(errfunc, p0[:], args=(x_o_l_2, y_o_l_2))
print("Optimized params:", p1)
plt.plot(x_o_l_2, y_o_l_2, "ro", x_o_l_2, fitfunc(p1,x_o_l_2), "r-", c='blue', alpha=0.5)
plt.xlim(1,3)
plt.show()
Потенциалом Морзе зависимость энергии молекулы от длины связи СС описывается гораздо лучше